Ab initio investigation of the melting line of nitrogen at high pressure 
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Understanding the behavior of molecular systems under pressure is a fundamental problem in 
condensed matter physics. In the case of nitrogen, the determination of the phase diagram and 
in particular of the melting line, are largely open problems. Two independent experiments have 
reported the presence of a maximum in the nitrogen melting curve, below 90 GPa, however the 
position and the interpretation of the origin of such maximum differ. By means of ab initio molecular 
dynamics simulations based on density functional theory and thermodynamic integration techniques, 
we have determined the phase diagram of nitrogen in the range between 20 and 100 GPa. We find 
a maximum in the melting line, related to a transformation in the liquid, from molecular N2 to 
polymeric nitrogen accompanied by an insulator-to-metal transition. 

PACS numbers: 



The properties of elemental materials under pres- 
sure attract considerable attention in condensed matter 
physics, geophysics and planetary science^ 2 -. In partic- 
ular, nitrogen, with its intricate phase diagram, and its 
potential applications as an energetic material, has been 
widely studied in recent years; however its phase dia- 
gram, including its melting line as a function of pressure, 
is still the subject of heated debate. 

In its solid form nitrogen remains molecular up to rel- 
atively high pressure (P ~ 100 GPa) and its phase dia- 
gram exhibits a variety of competing crystalline phases 3 . 
At low P, N2 molecules interact via weak dipolar and 
quadrupolar forces, while N atoms are held together by 
a triple bond, which is the strongest chemical bond in 
nature. The ability of transforming the triple bonds 
of molecular nitrogen into single bonds would open the 
way to storing energy at very high density — . This is 
in principle possible: by pressurizing nitrogen to about 
110 GPa, non- molecular crystalline and/or amorphous 
phases are formed 3 ^, as predicted by pioneering theoret- 
ical works^ 7 -. A crystalline covalent form, dubbed "cubic 
gauche" and proposed theoretically 8 , was obtained and 
fully characterized^ in diamond anvil cell experiments. 
In this insulating phase every N atom forms three sin- 
gle covalent bonds with its neighbors arranged in a cubic 
lattice. Further covalent crystalline forms have been pre- 
dicted to occur at even higher P 10 i n , but are yet to be 



found in experiments. 

By analogy with high pressure solid phases, the exis- 
tence of non-molecular, liquid nitrogen was suggested as 
welU^, and very recently a first order liquid-liquid phase 
transition has been proposed^ 3 -, between a low density 
liquid molecular phase (LDL) and a high density liquid 
polymeric phase (HDL), located between 2000 and 6000 
K at ~ 80 GPa. Such structural transformation is accom- 
panied by metallization of fluid nitrogen, as observed in 
shock reverberation experiments^. Though uncommon 
in elemental liquids, a first-order liquid-liquid (LL) phase 
transition, from a molecular to an atomic phase, has been 
observed in phosphorus^, which is isovalent to nitrogen. 
Support in favor of a LL phase transition in nitrogen 
comes from the observation of a maximum in the melt- 
ing curv o 16 ' 17 , whose presence may be an indication of a 
change in the liquid properties 15 i 18 i 19 . A negative slope 
of the melting line may also be associated with the pres- 
ence of open crystalline structures, as e.g. in carbon and 
water, or with changes in the electronic structure of the 
system, for example a metal-insulator transition^, or pro- 
motions of valence electrons to electronic orbitals higher 
in energies than those occupied at low P, as found, e.g. 
in Cesium. 

The position and the character of the maximum in 
the melting curve of nitrogen are still matter of de- 
bat o 16 i 17 ' 20 i 21 According to Ref ^ the maximum is sharp 



and located at 50 GPa, and may possibly be the sig- 
nature of a triple point associated to a first order LL 
phase transition. Goncharov et a/ . 17 ' 20 measured instead 
a slight change in the melting line slope near 70 GPa. In 
addition, by performing in situ Raman scattering, they 
found no evidence of a LL phase transition, and related 
the maximum in the melting curve to polymorphic tran- 
sitions between crystalline molecular phases. 

In this Letter we report the theoretical melting line 
of nitrogen between 20 and 100 GPa as obtained from 
first principle molecular dynamics (MD) simulations. We 
predict that the melting temperature reaches a maximum 
between 80 and 90 GPa, in correspondence to a transition 
in the liquid phase involving both a structural modifica- 
tion from a molecular to a polymeric fluid, and a semi- 
conductor to metal transition. We show that close to the 
maximum, the liquid polymerizes and becomes denser 
than the corresponding molecular solid, thus giving rise 
to a negative slope in the PT melting curve. 

Calculations of melting lines can be obtained either 
by the two-phase simulation metho d 18 ' 22 ' 23 or by ther- 
modynamic integration (TI). The two approaches are in 
principle equivalent^!, but the two-phase method may 
require larger simulation cells and longer runs to achieve 
accuracy comparable to TI. We therefore determined the 
melting temperature of nitrogen at several different P 
by computing free energy differences between liquid and 
crystalline phases by TI, in a manner similar to previ- 
ous studiesi of C and Si melting line o 25 ' 26 . Our compu- 
tational framework relies upon Born-Oppenheimer MD 
simulation :) 27 ' 28 , where the electronic structure is solved 
within density functional theory (DFT). We used a gener- 
alized gradient approximation, PBE 29 , for the exchange 
and correlation functional, norm conserving pseudopo- 
tentials and a plane- wave expansion of the electronic or- 
bitals with a kinetic-energy cut-off of 60 Ry. We sim- 
ulated nitrogen in supercells containing 128 atoms with 
T-point sampling of the supercell Brillouin zone. The 
MD equations of motion are integrated with a time-step 
of 20 a.u., the temperature is controlled by stochastic 
velocity rescaling 3 -^ and the pressure is kept constant by 
first-order cell dynamics. This scheme yields good agree- 
ment with the P(V) curve of liquid N at 2000 K in Ref. 13 . 

The use of TI requires the availability of a potential 
to describe a reference system, for which we have chosen 
a classical force-field that provides a good description of 
molecular nitrogen at low Pi 3138 . Our TI protocol con- 
sists of three steps: (i) We compute the melting temper- 
ature of the reference system (T^ ) at a given pressure 
P by a two-phase simulation, using a supercell with 3400 
molecules and a simulation time of 100 ps. We note that 
the melting line of the classical system is monotonic as 
a function of P and does not exhibit any maximum, (ii) 
The (Gibbs) free energy differences (AG) of the refer- 
ence and the DFT/PBE systems are computed both for 
the solid and the liquid phase by adiabatic switch 3 - 2 - (AS) 
in 2 ps runs. The convergence of AG as a function of 
switching time has been tested by performing longer AS 
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FIG. 1: (color online) Proposed phase diagram of nitro- 
gen. The computed points on the melting line are indicated 
with blue circles and the slopes obtained from the Clausius- 
Clapeyron equation with blue straight lines. The experimen- 
tal points from Goncharov et a/.— are shown as solid dia- 
monds and those measured by Mukherjee and Boehlpi^ as 
empty diamonds. 



runs at 40 GPa and error bars on AG are obtained as 
the standard deviation of four independent runs. Longer 
AS runs, up to 10 ps have been performed to obtain an 
accurate evaluation of the free energy of the liquid at 90 
GPa, given the inability of the classical reference model 
to describe the dissociated liquid. The free energy dif- 
ferences between solid and liquid at different volumes 
need to be corrected because of the finite plane wave 
cutoff used in our simulations 25 . A correction term is 
computed by performing single point calculations for se- 
lected MD snapshots with a plane wave cutoff of 160 Ry, 
that is sufficient to yield well converged values of P at all 
volumes considered here, (iii) The melting temperature 
(T m ) of nitrogen at the PBE level is finally computed by 
reversible scalin g 33 ' 3439 : T m is located at the intersec- 
tion of the Gibbs free energy curves (G(T)) of the liquid 
and the solid phases, with an initial offset determined by 
AQ at T^* , obtained in step (ii). In addition, we have 
computed the slope of the melting curve at 65, 80.5 and 
90 GPa using the Clausius-Clapeyron equation^ 5 ., eval- 
uating the differences in density and enthalpy over 5 ps 
long independent ab initio MD runs, carried out in the 
canonical ensemble at constant pressure (NPT). 

Using the procedure discussed above, we have com- 
puted the melting temperature at five different values of 
P: the results are reported in Fig. [1] and compared with 
the experiments from Refsi 16 ' 17 . Our results agree within 
one error bar with the melting points measured by Gon- 
charov et aZiii, while they are not compatible with the 
presence of a cusp at 50 GPa, as in Ref J^. Our computed 
melting curve displays a maximum between 80.5 and 90 
GPa; the presence of a maximum is further supported by 
the opposite signs of the slope of the curve, computed 
at 80.5 and 90 GPa. The position of the maximum in 
the melting line is shifted toward higher pressure with 



respect to the measurements of RefUi by about 10 GPa. 

Our simulations show that the thermodynamically sta- 
ble crystalline phase between 80 and 90 GPa is molecular, 
even at T close to the melting point, in agreement with 
the experimental observations in iLi£. Yet, at variance 
with RefJi, we could not locate a polymorphic phase 
transition between the S and the e phase within this pres- 
sure range. By direct ah initio MD simulations in the 
NPT ensemble we observe that at 1500 K and 90 GPa 
phase e transform rapidly into 5 and we estimate the 
5/e phase boundary at P above 120 GPa at 1500 K. Ex- 
periments confirm that solid nitrogen is molecular in the 
(P,T) range of interest (80-90 GPa, - 1750K), however 
the presence of several competing metastable structures 
and strong hysteresis effects make the determination of 
the stable crystalline phase uncertain 3 - 6 -. Nevertheless the 
density differences between the various molecular solids 
observed at this P are small, compared to the density 
difference between the liquid and the solid phase; thus 
uncertainties in the relative stability of the molecular 
polymorphs may only result in small quantitative vari- 
ations of the melting line predicted by our simulations. 

Therefore we conclude that the presence of a maximum 
in the melting line stems from structural and electronic 
transformations occurring in the liquid, rather than in 
the solid phase. As first observed in RefJ^, and con- 
firmed by our simulations, liquid nitrogen undergoes a 
transition from a molecular to a polymeric phase, which 
at 2000 K occurs at ^88 GPa. The analysis of our liquid 
sample at 90 GPa shows the coexistence of molecular N2 
and chains of N atoms where triple bonds give way to 
longer single bonds, whose signature appears as a second 
peak in the radial distribution function (not shown) at 
^1.3 A, while the triple bond yields a peak at 1.1 A. The 
distribution of the bond angles around the tetrahcdral 
value (109.3) and the observed tendency to form pen- 
tagonal rings are signatures of sp 3 hybridization of the 
atomic orbitals, analogous to the one observed in cova- 
lent crystalline phases^^ of nitrogen. The formation of 
covalent chains causes a drop in the volume of the liquid, 
which becomes denser than its crystalline counterpart. 
At 90 GPa and at the predicted melting temperature of 
1705 K, liquid N is 1% denser than the solid. Such den- 
sity difference results in a negative slope of the melting 
line. 

If the molecular to polymeric LL transition was first 
order, as suggested in ReL^ 3 -, then the maximum in the 
melting curve predicted by our simulations would coin- 
cide with a triple point, and it would be a cusp (i.e. the 
melting line would have discontinuous derivatives at the 
maximum), as Ref. 16 . However the location of the maxi- 
mum found here is different from the one found in Ref. 16 
(at about 50 GPa) . If the LL transition was instead sec- 
ond order, then the derivative of the melting line would 
be a continuous function, similar to what observed, for 
example, in liquid carbo n 22 ! 26 . The two possible scenar- 
ios are illustrated in Fig. [3] We note that the charac- 
terization of the LL phase transition as first-order— is 
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FIG. 2: (color online). Computed melting line of nitrogen 
(see Fig.l) and liquid- liquid phase boundary. The red solid 
line (guide to the eye) indicates the presence of a triple point, 
occurring in the case of a first-order LL phase transition. In- 
stead, the green dotted line (guide to the eye) does not show 
any triple point, corresponding to a second order LDL-HDL 
transition. 
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FIG. 3: (color online) Computed electronic density of states 
for the liquid phase at 80, 90 and 120 GPa and for the Cubic 
Gauche (CG) crystalline phase. 



plausible but not definitive: the use of small periodic sim- 
ulation cells may make a second-order phase transition 
resemble a first-order one, and a statistically significant 
size scaling analysis was not performed^ 3 -. Unfortunately 
such analysis would likely require cells with at least 500 
molecules and, most importantly, simulation times of the 
order of several ns, that are outside the reach of current 
ab initio simulation techniques. 

As a consequence of the structural changes occurring 
upon compression, the electronic structure of the liquid 
undergoes major modifications (Fig. [3]). Up to 80 GPa 
liquid nitrogen is an insulator with a DFT/PBE gap of 
3.1 eV. At 90 GPa, the formation of chains and pentago- 
nal rings give rise to the appearance of defect-like states 
in the middle of the electronic gap. These states have 
anti-bonding n* character and are delocalized over the 




FIG. 4: (color online) Isosurface of the square modulus of the 
wavefunction of a single particle state with energy inside the 
electronic gap of liquid nitrogen at 90 GPa. 



polymeric chains (Fig. 0]). As P is further increased, 
the liquid loses its molecular character, an increasing 
number of polymeric chains are formed, and eventually 
metallization occurs. The electronic density of states of 
liquid nitrogen at 120 GPa shows indeed no gap. The 
observed metallization is consistent with an increase in 
electrical conductivity measured in shock reverberation 
experiments 14 . The density of states of the CG phase 
at 120 GPa is shown in the upper panel of Fig. [3] it 
is remarkable that the stable crystalline covalent poly- 



morphs of nitrogen are semiconducting (or insulating) 
up to a pressure as high as 240 GPa 3 . However a chain- 
like metallic crystalline polymorph was predicted to have 
an enthalpy close to that of the CG phased. 

In summary, by means of ab initio MD simulations 
and thermodynamic integration we have determined the 
melting line (T m (P)) of nitrogen up to 90 GPa. We have 
found that T m (P) exhibits a maximum between 80 and 
90 GPa which is related to a structural transformation 
in the liquid, from a molecular to a polymeric phase. 
This transformation is accompanied by an insulator to 
metal transition. Our computed melting temperatures 
are in fair agreement with those determined in recent di- 
amond anvil cell experiments^, and not compatible with 
the data of RefJ^, where melting was established by vi- 
sual inspection. If the transformation observed in the 
liquid corresponds to a first order phase transition, as 
suggested in RefJ^, then the maximum found here will 
coincide with a triple point and thus a cusp in the melt- 
ing line, as proposed by Ref.— . However, the shift in the 
transition of the electronic properties of the liquid, which 
undergoes metallization at high pressure, with respect to 
its structural transition is an indication against a first or- 
der transition in the liquid phase. Work is in progress to 
compute spectroscopic observables capable of unequivo- 
cally identifying different liquid phases. 
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